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We consider a system of coupled classical harmonic oscillators with spatially fluctuating nearest- 
neighbor force constants on a simple cubic lattice. The model is solved both by numerically diag- 
onalizing the Hamiltonian and by applying the single-bond coherent potential approximation. The 
results for the density of states g(to) are in excellent agreement with each other. As the degree of 
ON ■ disorder is increased the system becomes unstable due to the presence of negative force constants. 

If the system is near the borderline of stability a low-frequency peak appears in the reduced density 
of states g(jjj)/u> as a precursor of the instability. We argue that this peak is the analogon of the 
"boson peak", observed in structural glasses. By means of the level distance statistics we show that 
the peak is not associated with localized states. 
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A ubiquitous and rather intriguing feature in the physics of glasses is the anomalous behavior of the low-frequency 
part of the vibration spectrum and the corresponding thermal properties [[O-pJ. While the origin of the linear low- 
tempature specific heat is commonly attributed to the existence of double- well potentials or two-level systems [Q, 
' there is still considerable debate about the so-called "boson peak" . This peak shows up in the density of states (DOS) 
. as an excess contribution, compared to the usual Debye behaviour ( g(u>) cx uj 2 ). It is taken as being responsible 

for a number of features observed in specific heat and thermal conductivity measurements, as well as in Raman or 
neutron scattering data 0-|^J§^l- The boson peak seems also to persist at elevated temperatures. Here its relation 
to the liquid-glass transition and the corresponding relaxation dynamics is a matter of discussion |p|-|l3||. Due to the 
development of new experimental techniques allowing to perform Brillouin scattering measurements in the THz range 
H |l7f , as well as pertinent molecular dynamics simulations (l8| ||^] the question concerning the nature of the modes 
^j" 1 , in the boson peak region has gained much additional interest recently. 

Several models have been formulated to explain the physical origin of the boson peak. In the soft-potential model 
p3| , ^4f the existence of anharmonic localized vibrations is postulated to be an intrinsic and essential property of 
amorphous systems. The scattering of propagating phonons by these modes is then taken as the origin for the excess 
DOS in the boson peak region. 

Another, almost orthogonal approach considers harmonic degrees of freedom solely. Here the scattering of phonons 
is associated with the disorder in the force constants, or, equivalently, spatially fluctuating elastic constants p5[ . 
The first contribution in this direction has been the phonon-fracton approach pq,Effl, in which harmonic vibrational 
excitations on a percolating lattice were considered. However, a numerical simulation of this model (2^] showed that 
the crossover from propagating modes (phonons) to localized fractal excitations (fractons) does not lead to an excess 
DOS, although calculations using the coherent potential approximation (CPA) |2^,^j| had predicted such an excess. 
Moreover, there is no experimental evidence for spatially fractal behaviour of bulk glasses. 

More realistic disordered coupled harmonic oscillator models work in terms of a continuous distribution of force 
constants. Calculations of the vibration spectrum of such a model with a specific force constant distribution have been 
performed in effective medium approximation (EMA, amorphous version of the CPA) pfj|-p3[. In these calculations 
an excess contribution in the DOS appears in the frequency regime, where the phonon mean free path becomes 
comparable to the wavelength. However, for this model the existence of the excess DOS has not yet been checked 
numerically. 

Another interesting contribution to the boson peak discussion p"2| comes from the mode-coupling theory |34| of the 
liquid-glass transition. In the idealized glass state a maximum in the density fluctuation spectrum at finite frequencies 
is predicted which is associated with inhomogeneously broadened harmonic oscillations. Experimental measurements 
exhibiting a crossover from relaxational to vibrational dynamics || have been described successfully in terms of this 
theory |L3|] 

Calculations for harmonic vibrational excitations in liquids, performed numerically ^||,^9), as well in EMA [j36| , 
can hardly contribute to the present debate, because the Hamiltonians considered are instable [EtJ 
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In the present letter we investigate the properties of a simple model for vibrational excitations of disordered systems 
and show that excess low-frequency modes are present if the system is almost unstable. The model consists of a set of 
coupled scalar harmonic oscillators placed on a three-dimensional simple cubic lattice with lattice constant a = l. The 
coupling among the oscillators is modelled by nearest-neighbor force constants Kij , which are treated as independent 
(quenched) random variables, chosen according to a distribution with density P(_Ky). The corresponding Hamiltonian 
has off-diagonal elements Hij = and diagonal elements Tin = — Kij ■ For a stable system all eigenvalues 
Ai = — ojf are negative, cj, denoting the vibrational eigenfrequencies. In all calculations we have chosen a truncated 
Gaussian distribution with P{K) = Poexp{— (K — K ) 2 /2a 2 }6(K — K m { n ). Here 9{x) denotes the step function, P 
is a normalization constant, K and er denote the maximum value and the width, respectively. In our calculations 
we set K = 1. The lower cut-off of the force constant distribution is denoted by K m [ n . This cut-off is introduced to 
allow for the study of strongly disordered systems with a restricted amount of negative force constants. 

The model defined in this way has been analyzed, both, numerically by diagonalizing the Hamiltonian matrix, as 
well as in CPA. The DOS obtained from both methods are in excellent agreement with each other, except near the 
point of instability. 

For the numerical treatment we considered a cubic box of size L and imposed periodic boundary conditions. The 
resulting L 3 x L 3 Hamiltonian matrix was diagonalized using the NAG-LAPACK routine DSYEV. The distribution 
of eigenvalues exhibits gaps due to finite-size effects. In order to eliminate these effects we calculated the integrated 
density of levels Fl(X) = J2i ^(^ — ^i) f° r a given size L, smoothed this function and averaged it over different sizes, 
ranging from L = 10 to L = 14. This procedure yielded a function F(X), the derivative of which gives the density of 
levels n(A) = dF(A)/dA. From this the vibrational DOS follows as g(uj) — 2ujn(—uj 2 ). 

For an approximate solution of our model we have used the single-bond CPA p8|-f40||. This theory is for- 
mulated in terms of a frequency-dependent force constant ("self energy") T(u>), which can be visualized as the 
inverse of an acoustic dielectric function. The frequency-dependent complex sound velocity v(lo) is given by 
T(u>) — v(ui) 2 . In CPA the quantity T ist determined self-consistently in terms of the local Green's function 
Gq(z) = X) \ z + ^ — 2(cos k x + cos k y + cos fc z )] _1 of the ordered cubic lattice (the sum runs over the Brillouin 

zone and z = —uj 2 + ie) as follows: 



T(z) - Kij 



1 - (r(z) - Kj) (l - zGo(z)) /3T(z) , 



(1) 



with G (z) = G (z/T(z))/T(z) and 

(A) = J &KijP{Kij)A(Kij). The DOS is obtained from a numerical solution of eq. ([!]) and use of the relation 
g(w) = -{2/n)w1m{G(z)}. 

Preliminary calculations showed that without lower cutoff (i. e. K mm = — oo) the system becomes instable for 
a w 0.6. In order to be able to study a strongly disordered system we have set a = 1 and varied the cutoff K m i n 
which controls the amount of small and negative force constants. Fig. 1 shows results of such calculations. The 
excellent agreement of the CPA calculations with the numerical analysis [[llj (except for the immediate vicinity of the 
instability) indicates both the reliability of the CPA and the correctness of the procedure utilized for the elimination 
of finite-size effects in the numerical work. As K m ^ n decreases from positive to negative values, the peak in g(u))/u) 2 
shifts towards smaller frequencies. At the same time the peak intensity increases. In CPA the system becomes 
instable for K m [ n = —0.85, whereas in the numerical calculation the instability already occurs near K m [ n = —0.6. 
The low-frequency peak thus plays the role of a precursor of the instability introduced by the presence of small 
and negative force constants. We therefore conclude that a low-frequency peak observed in the reduced DOS of a 
disordered harmonic system indicates that it is almost unstable. 

Another quantity of interest, in particular in the context of thermal conductivity data in structural glasses [p|,^2[, 
is the mean free path of the phonons I. Within the CPA treatment the mean free path can be identified with the 
decay length of the wave intensity | exp{iu>r /v(u>)}\ 2 and is therefore given by 



1 {UJ) ~ \v(u)\ 2 [) 
and behaves as l~ x oc w 4 for lu — > (Rayleigh scattering). 

In fig. 2 we have plotted the reduced quantity € _1 (cu) / 'uo A for the same force constant distributions as in fig. 1. It is 
seen that the mean free path strongly decreases with the amount of small and negative force constants. In scattering 
experiments the decrease in l(uS) would show up as an increasing width of the Brillouin peaks. The uj a dependence of 
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£(u}) levels off near the frequency where the reduced DOS has its maximum. It has already been demonstrated in the 
literature that such a frequency dependence of £ is capable to explain the plateau in the thermal conductivity ^,^2) . 
One sees from fig. 2 that for K m { n = —0.6 (almost unstable case) I is in this frequency range of order unity. Since 
the real part of v{u) and uj are of order unity as well the mean free path becomes comparable to the wavelength, in 
which case the phonons may become localized according to the Ioffe-Regel criterion p3| . 

Therefore we were interested in the question as to whether in our model the vibrational states are localized. This 
can be checked by an investigation of the level statistics p3| . If we define normalized eigenvalues by = F{\i) a 
histogram of the distances Si — |ej — should yield a distribution according to the Gaussian orthogonal random 
matrix ensemble (GOE) P(s) = ^irs exp{— 7rs 2 /4} in the case of delocalized states, whereas one expects a Poissonian 
distribution P(s) — exp{ — s} for localized states. This is due to the fact, that the delocalized states show level 
repulsion, whereas the localized ones do not. 

In fig. 3 we show a representative set of the level distance statistics of our calculated spectra. It is clearly seen that 
they universally follow the GOE statistics. We therefore conclude that the states in our model are all delocalized for 
the force constant distributions considered, also in the vicinity of the instability. 

When comparing our results with real three-dimensional disordered systems like glasses a few comments are in 
order. We have treated a scalar model and thus ignored the vector character of the phonons completely. We believe, 
however, that our model already is able to mimic the behavior of transverse phonons in a glass. If the vector character 
is taken into account, one would expect an even stronger scattering due to the admixture of longitudinal degrees of 
freedom. We therefore believe that our model underestimates the intensity of the boson peak somewhat. We expect 
also that in a modified version of our model the states near the boson peak may become Anderson-localized near 
the borderline of stability. Therefore we have to state that in a real glassy system it still remains an open question, 
whether the modes associated with the boson peak are localized or delocalized. Our model shows, however, that is not 
necessary to postulate the existence of localized states or strongly anharmonic effects to obtain an excess contribution 
in the DOS. 

As we have treated only harmonic interactions the effects associated with the glass-liquid transformation [jl2 34 
are outside the scope of the present model. 

In conclusion we have solved a simple model of coupled harmonic oscillators, both, numerically and in coherent- 
potential approximation. Near the point of instability the model exhibits a low-frequency peak in the reduced density 
of states g{uj)/w 2 , which we view as the analogue to the boson peak observed in structural glasses. 

It is a pleasure to thank W. Gotze, W. Petry, A. Meyer, and J. Wuttke for fruitful discussions and a critical reading 
of the manuscript. 
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FIG. 1: 

Reduced DOS g(to)/u} 2 against frequency for force constant distributions with a = Kq = 1 and different lower cutoffs 
if m i n . The symbols represent the numerical diagonalization, the full lines the CPA results. The agreement is achieved 
without any adjustable parameters. 

FIG. 2: 

Reduced inverse mean free path l~ x /uj 4 calculated in CPA according to (j^) for the parameters of fig. 1. 
FIG. 3: 

Level distance statistics for the distributions with negative values of K m [ n . Since they follow the statistics of the 
Gaussian orthogonal ensemble (GOE), the corresponding states are delocalized. 
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Fig. 2 
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